/********************************************************************
 *                                                                  *
 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
 *                                                                  *
 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2002             *
 * by the XIPHOPHORUS Company http://www.xiph.org/                  *
 *                                                                  *
 ********************************************************************

 function: PCM data envelope analysis
 last mod: $Id: envelope.c,v 1.54 2003/09/05 23:17:49 giles Exp $

 ********************************************************************/

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "../ogg/ogg.h"
#include "../vorbis/codec.h"
#include "codec_internal.h"

#include "os.h"
#include "scales.h"
#include "envelope.h"
#include "mdct.h"
#include "misc.h"

void _ve_envelope_init( envelope_lookup* e, vorbis_info* vi )
{
	codec_setup_info* ci = vi->codec_setup;
	vorbis_info_psy_global* gi = &ci->psy_g_param;
	int ch = vi->channels;
	int i, j;
	int n = e->winlength = 128;
	e->searchstep = 64; /* not random */

	e->minenergy = gi->preecho_minenergy;
	e->ch = ch;
	e->storage = 128;
	e->cursor = ci->blocksizes[1] / 2;
	e->mdct_win = _ogg_calloc( n, sizeof( *e->mdct_win ) );
	mdct_init( &e->mdct, n );

	for( i = 0; i < n; i++ )
	{
		e->mdct_win[i] = sin( i / ( n - 1. ) * M_PI );
		e->mdct_win[i] *= e->mdct_win[i];
	}

	/* magic follows */
	e->band[0].begin = 2;
	e->band[0].end = 4;
	e->band[1].begin = 4;
	e->band[1].end = 5;
	e->band[2].begin = 6;
	e->band[2].end = 6;
	e->band[3].begin = 9;
	e->band[3].end = 8;
	e->band[4].begin = 13;
	e->band[4].end = 8;
	e->band[5].begin = 17;
	e->band[5].end = 8;
	e->band[6].begin = 22;
	e->band[6].end = 8;

	for( j = 0; j < VE_BANDS; j++ )
	{
		n = e->band[j].end;
		e->band[j].window = _ogg_malloc( n * sizeof( *e->band[0].window ) );
		for( i = 0; i < n; i++ )
		{
			e->band[j].window[i] = sin( ( i + .5 ) / n * M_PI );
			e->band[j].total += e->band[j].window[i];
		}
		e->band[j].total = 1. / e->band[j].total;
	}

	e->filter = _ogg_calloc( VE_BANDS * ch, sizeof( *e->filter ) );
	e->mark = _ogg_calloc( e->storage, sizeof( *e->mark ) );

}

void _ve_envelope_clear( envelope_lookup* e )
{
	int i;
	mdct_clear( &e->mdct );
	for( i = 0; i < VE_BANDS; i++ )
	{
		_ogg_free( e->band[i].window );
	}
	_ogg_free( e->mdct_win );
	_ogg_free( e->filter );
	_ogg_free( e->mark );
	memset( e, 0, sizeof( *e ) );
}

/* fairly straight threshhold-by-band based until we find something
   that works better and isn't patented. */

static int _ve_amp( envelope_lookup* ve,
					vorbis_info_psy_global* gi,
					float* data,
					envelope_band* bands,
					envelope_filter_state* filters,
					long pos )
{
	long n = ve->winlength;
	int ret = 0;
	long i, j;
	float decay;

	/* we want to have a 'minimum bar' for energy, else we're just
	   basing blocks on quantization noise that outweighs the signal
	   itself (for low power signals) */

	float minV = ve->minenergy;
	float* vec = alloca( n * sizeof( *vec ) );

	/* stretch is used to gradually lengthen the number of windows
	   considered prevoius-to-potential-trigger */
	int stretch = max( VE_MINSTRETCH, ve->stretch / 2 );
	float penalty = gi->stretch_penalty - ( ve->stretch / 2 - VE_MINSTRETCH );
	if( penalty < 0.f )
	{
		penalty = 0.f;
	}
	if( penalty > gi->stretch_penalty )
	{
		penalty = gi->stretch_penalty;
	}

	/*_analysis_output_always("lpcm",seq2,data,n,0,0,
	  totalshift+pos*ve->searchstep);*/

	/* window and transform */
	for( i = 0; i < n; i++ )
	{
		vec[i] = data[i] * ve->mdct_win[i];
	}
	mdct_forward( &ve->mdct, vec, vec );

	/*_analysis_output_always("mdct",seq2,vec,n/2,0,1,0); */

	/* near-DC spreading function; this has nothing to do with
	   psychoacoustics, just sidelobe leakage and window size */
	{
		float temp = vec[0] * vec[0] + .7 * vec[1] * vec[1] + .2 * vec[2] * vec[2];
		int ptr = filters->nearptr;

		/* the accumulation is regularly refreshed from scratch to avoid
		   floating point creep */
		if( ptr == 0 )
		{
			decay = filters->nearDC_acc = filters->nearDC_partialacc + temp;
			filters->nearDC_partialacc = temp;
		}
		else
		{
			decay = filters->nearDC_acc += temp;
			filters->nearDC_partialacc += temp;
		}
		filters->nearDC_acc -= filters->nearDC[ptr];
		filters->nearDC[ptr] = temp;

		decay *= ( 1. / ( VE_NEARDC + 1 ) );
		filters->nearptr++;
		if( filters->nearptr >= VE_NEARDC )
		{
			filters->nearptr = 0;
		}
		decay = todB( &decay ) * .5 - 15.f;
	}

	/* perform spreading and limiting, also smooth the spectrum.  yes,
	   the MDCT results in all real coefficients, but it still *behaves*
	   like real/imaginary pairs */
	for( i = 0; i < n / 2; i += 2 )
	{
		float val = vec[i] * vec[i] + vec[i + 1] * vec[i + 1];
		val = todB( &val ) * .5f;
		if( val < decay )
		{
			val = decay;
		}
		if( val < minV )
		{
			val = minV;
		}
		vec[i >> 1] = val;
		decay -= 8.;
	}

	/*_analysis_output_always("spread",seq2++,vec,n/4,0,0,0);*/

	/* perform preecho/postecho triggering by band */
	for( j = 0; j < VE_BANDS; j++ )
	{
		float acc = 0.;
		float valmax, valmin;

		/* accumulate amplitude */
		for( i = 0; i < bands[j].end; i++ )
		{
			acc += vec[i + bands[j].begin] * bands[j].window[i];
		}

		acc *= bands[j].total;

		/* convert amplitude to delta */
		{
			int p, this = filters[j].ampptr;
			float postmax, postmin, premax = -99999.f, premin = 99999.f;

			p = this;
			p--;
			if( p < 0 )
			{
				p += VE_AMP;
			}
			postmax = max( acc, filters[j].ampbuf[p] );
			postmin = min( acc, filters[j].ampbuf[p] );

			for( i = 0; i < stretch; i++ )
			{
				p--;
				if( p < 0 )
				{
					p += VE_AMP;
				}
				premax = max( premax, filters[j].ampbuf[p] );
				premin = min( premin, filters[j].ampbuf[p] );
			}

			valmin = postmin - premin;
			valmax = postmax - premax;

			/*filters[j].markers[pos]=valmax;*/
			filters[j].ampbuf[this] = acc;
			filters[j].ampptr++;
			if( filters[j].ampptr >= VE_AMP )
			{
				filters[j].ampptr = 0;
			}
		}

		/* look at min/max, decide trigger */
		if( valmax > gi->preecho_thresh[j] + penalty )
		{
			ret |= 1;
			ret |= 4;
		}
		if( valmin < gi->postecho_thresh[j] - penalty )
		{
			ret |= 2;
		}
	}

	return( ret );
}

#if 0
	static int seq = 0;
	static ogg_int64_t totalshift = -1024;
#endif

long _ve_envelope_search( vorbis_dsp_state* v )
{
	vorbis_info* vi = v->vi;
	codec_setup_info* ci = vi->codec_setup;
	vorbis_info_psy_global* gi = &ci->psy_g_param;
	envelope_lookup* ve = ( ( private_state* )( v->backend_state ) )->ve;
	long i, j;

	int first = ve->current / ve->searchstep;
	int last = v->pcm_current / ve->searchstep - VE_WIN;
	if( first < 0 )
	{
		first = 0;
	}

	/* make sure we have enough storage to match the PCM */
	if( last + VE_WIN + VE_POST > ve->storage )
	{
		ve->storage = last + VE_WIN + VE_POST; /* be sure */
		ve->mark = _ogg_realloc( ve->mark, ve->storage * sizeof( *ve->mark ) );
	}

	for( j = first; j < last; j++ )
	{
		int ret = 0;

		ve->stretch++;
		if( ve->stretch > VE_MAXSTRETCH * 2 )
		{
			ve->stretch = VE_MAXSTRETCH * 2;
		}

		for( i = 0; i < ve->ch; i++ )
		{
			float* pcm = v->pcm[i] + ve->searchstep * ( j );
			ret |= _ve_amp( ve, gi, pcm, ve->band, ve->filter + i * VE_BANDS, j );
		}

		ve->mark[j + VE_POST] = 0;
		if( ret & 1 )
		{
			ve->mark[j] = 1;
			ve->mark[j + 1] = 1;
		}

		if( ret & 2 )
		{
			ve->mark[j] = 1;
			if( j > 0 )
			{
				ve->mark[j - 1] = 1;
			}
		}

		if( ret & 4 )
		{
			ve->stretch = -1;
		}
	}

	ve->current = last * ve->searchstep;

	{
		long centerW = v->centerW;
		long testW =
			centerW +
			ci->blocksizes[v->W] / 4 +
			ci->blocksizes[1] / 2 +
			ci->blocksizes[0] / 4;

		j = ve->cursor;

		while( j < ve->current - ( ve->searchstep ) )
		{
			/* account for postecho
			                                         working back one window */
			if( j >= testW )
			{
				return( 1 );
			}

			ve->cursor = j;

			if( ve->mark[j / ve->searchstep] )
			{
				if( j > centerW )
				{

#if 0
					if( j > ve->curmark )
					{
						float* marker = alloca( v->pcm_current * sizeof( *marker ) );
						int l, m;
						memset( marker, 0, sizeof( *marker )*v->pcm_current );
						fprintf( stderr, "mark! seq=%d, cursor:%fs time:%fs\n",
								 seq,
								 ( totalshift + ve->cursor ) / 44100.,
								 ( totalshift + j ) / 44100. );
						_analysis_output_always( "pcmL", seq, v->pcm[0], v->pcm_current, 0, 0, totalshift );
						_analysis_output_always( "pcmR", seq, v->pcm[1], v->pcm_current, 0, 0, totalshift );

						_analysis_output_always( "markL", seq, v->pcm[0], j, 0, 0, totalshift );
						_analysis_output_always( "markR", seq, v->pcm[1], j, 0, 0, totalshift );

						for( m = 0; m < VE_BANDS; m++ )
						{
							char buf[80];
							sprintf( buf, "delL%d", m );
							for( l = 0; l < last; l++ )
							{
								marker[l * ve->searchstep] = ve->filter[m].markers[l] * .1;
							}
							_analysis_output_always( buf, seq, marker, v->pcm_current, 0, 0, totalshift );
						}

						for( m = 0; m < VE_BANDS; m++ )
						{
							char buf[80];
							sprintf( buf, "delR%d", m );
							for( l = 0; l < last; l++ )
							{
								marker[l * ve->searchstep] = ve->filter[m + VE_BANDS].markers[l] * .1;
							}
							_analysis_output_always( buf, seq, marker, v->pcm_current, 0, 0, totalshift );
						}

						for( l = 0; l < last; l++ )
						{
							marker[l * ve->searchstep] = ve->mark[l] * .4;
						}
						_analysis_output_always( "mark", seq, marker, v->pcm_current, 0, 0, totalshift );


						seq++;

					}
#endif

					ve->curmark = j;
					if( j >= testW )
					{
						return( 1 );
					}
					return( 0 );
				}
			}
			j += ve->searchstep;
		}
	}

	return( -1 );
}

int _ve_envelope_mark( vorbis_dsp_state* v )
{
	envelope_lookup* ve = ( ( private_state* )( v->backend_state ) )->ve;
	vorbis_info* vi = v->vi;
	codec_setup_info* ci = vi->codec_setup;
	long centerW = v->centerW;
	long beginW = centerW - ci->blocksizes[v->W] / 4;
	long endW = centerW + ci->blocksizes[v->W] / 4;
	if( v->W )
	{
		beginW -= ci->blocksizes[v->lW] / 4;
		endW += ci->blocksizes[v->nW] / 4;
	}
	else
	{
		beginW -= ci->blocksizes[0] / 4;
		endW += ci->blocksizes[0] / 4;
	}

	if( ve->curmark >= beginW && ve->curmark < endW )
	{
		return( 1 );
	}
	{
		long first = beginW / ve->searchstep;
		long last = endW / ve->searchstep;
		long i;
		for( i = first; i < last; i++ )
			if( ve->mark[i] )
			{
				return( 1 );
			}
	}
	return( 0 );
}

void _ve_envelope_shift( envelope_lookup* e, long shift )
{
	int smallsize = e->current / e->searchstep + VE_POST; /* adjust for placing marks
						     ahead of ve->current */
	int smallshift = shift / e->searchstep;

	memmove( e->mark, e->mark + smallshift, ( smallsize - smallshift )*sizeof( *e->mark ) );

#if 0
	for( i = 0; i < VE_BANDS * e->ch; i++ )
		memmove( e->filter[i].markers,
				 e->filter[i].markers + smallshift,
				 ( 1024 - smallshift )*sizeof( *( *e->filter ).markers ) );
	totalshift += shift;
#endif

	e->current -= shift;
	if( e->curmark >= 0 )
	{
		e->curmark -= shift;
	}
	e->cursor -= shift;
}






